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Structural optimization with a flutter constraint for a vehicle designed to fly in the transonic regime is a 
particularly difficult task. In this speed range, the flutter boundary is very sensitive to aerodynamic nonlinear- 
ities, typically requiring high-fidelity Navier-Stokes simulations. However, the repeated application of unsteady 
computational fluid dynamics to guide an aeroelastic optimization process is very computationally expensive. 
This expense has motivated the development of methods that incorporate aspects of the aerodynamic nonlin- 
earity, classical tools of flutter analysis, and more recent methods of optimization. While it is possible to use 
doublet lattice method aerodynamics, this paper focuses on the use of an unsteady high-fidelity aerodynamic 
reduced order model combined with successive transformations that allows for an economical way of utiliz- 
ing high-fidelity aerodynamics in the optimization process. This approach is applied to the common research 
model wing structural design. As might be expected, the high-fidelity aerodynamics produces a heavier wing 
than that optimized with doublet lattice aerodynamics. It is found that the optimized lower skin of the wing 
using high-fidelity aerodynamics differs significantly from that using doublet lattice aerodynamics. 


Nomenclature 
Variables 
coefficient matrix of dimension N x N, 


damping matrix 


Cc 

D 

F nodal force vector, N 

g generalized displacement vector, m 

G generalized force vector, N 

H(s) frequency domain aerodynamic influence coefficient (AIC) matrix of 
dimension N x N, m 

k reduced frequency, k = OLr/Uco 

K stiffness matrix, N/m 

Lr reference length 

M mass matrix, kg 

N number of node points 

N, number of retained modes 

Ne number of modes in the expanded set, N, = N, 

oo dynamic pressure, N / m 

O(s) frequency domain generalized aerodynamic influence coefficient (g AIC) 
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matrix of dimension N, x N,, m 


Dp nondimensional Laplace variable, p = sLr/Ueo 

s dimensional Laplace variable, s = Cw + ia, s ' 

Ug freestream velocity, m/s 

UeEas freestream equivalent air speed, m/s 

x vector of displacements, m 

n minimum modal flutter damping 

An n'" chosen root in the Roger approximation of the generalized AIC matrix 

o undamped natural frequency of structure, rad / sec 

Or frequency of flutter, rad / sec 

Q diagonalized matrix of free vibration frequencies of the structure, rad / sec 

C true damping, € = 1/22[/n(an+1/a,)], where a, is the n'" oscillation peak 
Superscripts 

~ Laplace transform 

* chosen optimization constants 
Subscripts 

b baseline 

m modified 


I. Introduction 

Aeroelastic optimization of flexible wing structures typically seeks to minimize structural weight and/or fuel burn, 
with constraints placed upon various flight metrics. Of these constraints, the aeroelastic flutter margin is particularly 
critical, due to the abrupt, damaging, and sensitive nature of this failure mechanism, and also because a large fraction of 
the structural weight is typically allocated to satisfying the flutter constraint [1]. For transonic applications, the flutter 
boundary is very sensitive to aerodynamic nonlinearities, typically requiring high-fidelity Navier-Stokes simulations 
for its accurate prediction [2]. At the same time, flutter simulations based on computational fluid dynamics (CFD) are 
very computationally expensive, and somewhat incompatible to the high volume of function evaluations required for 
structural and aeroelastic optimization in an unsteady flowfield. 

Aeroelastic optimization using Reynolds-averaged Navier-Stokes (RANS) CFD, nevertheless, is desirable since 
flowfield nonlinearity can have a significant effect on the optimized structure. With the progress in computing power, 
this approach to performing at least static aeroelastic optimization seems feasible for some cases, as attested to by 
recent published results [3—1 1]. In addition, there have been recent investigations into optimization of shape or aeroe- 
lastic properties using the adjoint method with the time-domain [1 2—15] or frequency-domain [16] Euler or Reynolds- 
averaged Navier-Stokes equations. In all of those time or frequency-domain studies, the configurations have either 
been two-dimensional or simple three-dimensional wings. 

Given the complexity of industrial flight configurations and the growing recognition of the amount of grid reso- 
lution required for solution accuracy, flutter- or gust-based optimization of a structure (or of a coupled aerodynamics, 
structure and control system) driven by the unsteady Reynolds-averaged Navier-Stokes (URANS) equations, still 
seems to be years away from routine use. This limitation encourages finding other means to incorporate high-fidelity 
aerodynamics into a structural optimization when constrained by a required flutter onset margin. 

Obviously, there has always been a motivation to understand and develop effective methods of including flutter 
constraints during structural optimization [17, 18]. Haftka and Yates [19] pioneered the development of a flutter 
constrained optimization method, using the low fidelity aerodynamic tools available at that time. More recent studies 


have addressed the problem of optimization with a flutter or gust load constraint [14, 20-26], but still used simplified 


aerodynamics such as Euler CFD, doublet lattice, unsteady panel, or strip theory aerodynamics. 

Linear unsteady aerodynamic tools exist that can compute an offline database of aerodynamic coefficients for 
a given planform [27—29]. As these coefficients are not functions of structural design variables, they may be used 
across the entire aeroelastic optimization process for an efficient handling of the flutter margin constraint. Linear 
flutter constraint boundaries are nonconservative (i.e., they overpredict the flutter dynamic pressure), and so a variety 
of tools exist to correct the unsteady aerodynamics, either based on CFD or wind tunnel data [30]. Overset-field panel 
methods have also been developed via oscillatory linearizations about a steady CFD result, interpolated onto a linear 
panel mesh [31]; the final result is a database of offline transonic aerodynamic coefficients, which as before, is highly 
amenable to aeroelastic optimization. Chen, et al. [32] have developed a method of transforming unsteady pressure 
coefficients from a baseline to a modified set of modes using a linearization of data from a nonlinear transonic small 
disturbance code. 

In this paper, a somewhat simpler and possibly less expensive formulation than that in Ref. 32 is developed to 
transform unsteady aerodynamic data to an updated structural model. This method provides an economical way to 
generate a CFD-based flutter margin. The unsteady data is derived by using aeroelastic system identification tools [33] 
to output a generalized aerodynamic influence coefficient (gAIC) matrix. The Roger approximation [34] of either 
the URANS or doublet lattice method (DLM) aerodynamics are combined with the equations of structural dynamics 
to obtain flutter predictions at select Mach numbers. Ostensibly, the Roger approximation is only useful for the 
structure whose mode shapes (eigenvectors) were used for structural perturbations to generate the CFD-based unsteady 
aerodynamics. Methods exist, however, to project the Roger approximation terms from the baseline modal space into 
a modified modal space, making the technique well-suited for aeroelastic optimization, where the mode shapes will 


change from one design iteration to the next. 


II. Aerodynamic Relationship between the Baseline and Modified Structures 

The present goal is to drive a structural model from a baseline configuration toward a modified optimal structure 
under a flutter constraint utilizing CFD-based unsteady aerodynamics. It has been reported that under certain condi- 
tions the influence of the derivatives of the mode shapes on the generalized aerodynamic force is negligible. In this 
case, the terms carrying the influence of mode shape change on the gAIC can be neglected and the optimization is 
influenced mostly by natural frequency of the optimized structure [20, 22]. In the present case where high-fidelity 
CFD at transonic flight conditions is used, it is unclear whether the effect of a mode shape change on the gAIC will 
be insignificant. Thus, a more general approach is likely to be needed. In the approach used here, the high-fidelity 
unsteady aerodynamics of the baseline model is used to derive the gAIC matrix of the modified structure. This is done 
via a projection of the modified modes onto the baseline modal basis vectors. 

Neglecting the structural damping, the equations of structural dynamics for the baseline structural model can be 


written 
M,X+K,px=F, . (1) 


The unknown x(r) is the vector of dynamic displacement increments and F(t) is the vector of aerodynamic forcing 


increments at N node points, both from a steady-state equilibrium. A Laplace transformation yields 
2 = ~ 
[s°Mp + Ky ]%(s) = Fy(s) (2) 


These equations are diagonalized with the mass normalized eigenvectors ®,, such that D) M,®, =I, ©) K,®, = QF 
and x = ®,g. Thus, 


[s°1+O5]8(s) =O, F,=Gy (3) 


There is a unique solution g, to this problem everywhere the left hand side is invertable. Note that functional depen- 
dence of g and Gon s will be implied but left out from here on. 

The optimization proceeds to a modified structural model with eigenvectors ®,,, such that O) My Pn =T7, 
©) Kn Pn = Qe and x = ®,,g and finally to an expression similar to equation 3. From here on, a unified notation 
will be used in which the subscript y is either b for baseline or m for modified model. With this notation, the aerody- 


namic forcing of the structure can be written 


Fy(s) = doo (s)3(s)y Y= (b,m) . (4) 


Since the AIC matrix H(s) is node dependent only, and absent large geometric changes, the same nodal AIC matrix 
will apply to the baseline and modified structures [19, 32]. 


Write equation 4 in generalized coordinates 
ss T T z e 
Gy = By Fy = qooPyH (5) PyBy = qooOy(5) By (5) 


where Qy is the gAIC matrix. A relationship between the gAIC of the baseline and modified structures is required. An 
orthogonal basis can be used such as that obtained by an application of the Gram-Schmidt procedure [35, 36] to the 
baseline modes. On the other hand, it is sometimes advantageous to use a nonorthogonal basis, for example, the use 
of linear normal modes as bases for a nonlinear structure [37-40]. In the present case, the modified structure mode 


shapes can be expressed in terms of the eigenvectors of the baseline undamped structure using 
@,=,C, . (6) 


Cy is a matrix of constant coefficients where C, = J for the baseline structure and C,, # J for the modified structure. 


See also equation 3 of Ref. 32. Substituting equation 6 into equation 5 yields 

~ T aT p T ss 

Gy = qooCy Pp HPCyBy = qooCy QpCyBy , Y= (bm), (7) 
which reveals that 


Ons) =CyOp(s)Cy (8) 


Equation 8 provides a means to calculate the gAIC of the modified structure from the baseline structure gAIC. 


To calculate C, Ref. 32 starts from equation 6 and uses a least-squares solution given by 
Ta \-l aT 
Cy = (Pp, Pp) OB, . (9) 
One can also compute the inverse of ®, if it is square, or pseudo-inverse otherwise, to get 
-1 
Cy=0;'0, . (10) 


An alternate approach is to use equation 6, the orthogonality of the baseline mode shapes with respect to the baseline 


mass matrix, and mass normalization of the baseline modal vectors to get 
Zr 
Cy =P,M,Py . (11) 


The most efficient method of computing C, can be determined on a case by case basis. 

When performing CFD-based flutter or system identification, the mode shapes of the structure at structural nodes 
must be projected onto the surface of the CFD model. In this context, it is important to note that in each case the 
modal displacement at the structural model nodes are used in equations 9- 11 rather than that projected to the CFD 
surface. Thus, the modified mode shapes do not need to be projected onto the CFD surface and likewise a repeated 


CFD analysis of the modified structure is not required. This simplification results in a significant reduction in the 
computation required. 

For a complete baseline eigenvector basis, the representation of the modified modes will be exact. For a model 
with N degrees of freedom, but with only N, modes retained (typically N, « N), the modified modes calculated by 
equation 6 will be approximate. To improve the accuracy requires the use of an enriched basis set to calculate C,,,. It is 
possible to get a reduced set of modes ®,,, that is exact by using all N basis vectors (®,) in calculating the coefficient 
matrix C,,, which in this case will have dimensions N x N,. From a practical standpoint, a few additional modes used 
to calculate Q,, can result in increased accuracy. One can partition the baseline gAIC Q, into an N,. x N, submatrix 
designated by subcript rr and additional submatrices re, er and ee (ee submatrix having dimensions (N, — N,) x (N. - 
N,.)) to get 


_ (Qp) rr (Qp)er 
o1=| [orn Hl | oy 


The addition of a few modes beyond N, may result in somewhat better accuracy in the retained modes without a large 
increase in the system identification of the N. mode responses. In this case C,, will have dimensions N x N,. 
A way to assess the adequacy of the NV, modes is to do a cross correlation of the approximated modified modes 


(o',) using equation 6 and the exact modified modes (®,,) 
I 
B = (On * Pn) / (Pn Pn) (13) 


where an exact representation will give B = 1. 


Il. Flutter Equation and Optimization Method 
The aerodynamic stiffness, damping, mass and wake induced lag effect are introduced by the Roger rational 
function approximation [34], which are transformed into the nondimensional Laplace variable p = sLr/Ueo 


Ni lags 


Ov(p)& = [Aon + PAiy + P Ary + >. PA(n+2)p/ (P — An) ]& (14) 


n=1 


where Njags is the number of lag states with root A, for lag n. These equations embody the unsteady aerodynamics 


which can be supplied by the DLM or URANS equations. The flutter equations (which now include a damping term) 


are written 
Weg 0 
co co oe 
La Map #5 Peto Te=0 (15) 
where 


MyaT Geol Ant, « DyS2bQysguC Ant, 
‘i ee (16) 
Ky = Qy = qooCy AopCy » F=oo y Cy PA(n+2),Cyl (P —A,)] 
n=l 

These equations are reformulated in state-space form, and the migration of the system’s Laplacian eigenvalues, 

as a function of dynamic pressure goo (or more commonly, equivalent air speed Ug,s), are simply computed with the 
p-method [41]. The Laplacian may be expanded as p = n +i-k; an 7 value greater than zero indicates an unstable 
system, and 7) is equal to zero at the flutter point. Rather than actually locate these flutter points, a method more 


amenable to the gradient based optimization used here is discussed by Hajela [42] and Ringertz [18], where the traces 


of 7 as a function of Ug4s must be below some bounding curve. This method allows the flutter mechanism to change 


identity during the optimization process, while still maintaining a smooth design space. Specifically: 


n*«(3-U* -Ugas—2+Usas)/(U*)? 0s Ugas < U* 


*\2 * * (17) 
c+ (Ugas -U~)* — Ugas > U 


where U~ is the minimum allowable flutter equivalent air speed, typically based on a 20% margin, 7” is the requisite 
damping at this flutter margin, and c is a quadratic scaling parameter. If 7* is set to 0, the method of Ringertz [18] 
is recovered. Critical Ug,s points (local minima) of the inequality in Eq. 17 are computed and lumped together into a 
single aggregate constraint via the Kreisselmeier-Steinhauser (KS) function [43]. 

For the gradient based optimization exercises discussed below, the derivative of critical points in Eq. 17 with 
respect to various structural sizing design variables must be computed analytically. The eigenvalue derivatives are 
simply computed using methods in Ref. 44 via a fixed-mode approximation, where the derivative of the mode shapes 
(®,,,) with respect to design variables is neglected for the purposes of gradient computations. 

During the optimization process (discussed in detail below), the numerical optimizer will continually update the 
structural details of the wing. At a given design iteration, the mode shapes of the current structure are designated ®,,,. 
The baseline mass matrix M,, mode shapes ®,, and gAIC matrix Q, are all constant and read from disk at the start of 
the procedure, and the procedures outlined above are used to compute the gAIC of the current structural iterate, Qy. 
The analytical derivatives of these baseline terms (Mz, ®,, Q,) with respect to the modified structural layout design 
variables are exactly zero. Derivatives of Qy are approximated to be zero via the fixed-mode approximation discussed 
above. 


IV. Results 

As an example, consider the undeflected common research model (uCRM). The uCRM is a jig shape of the lg 
wing shape model developed in Ref. 45. The uCRM is a generic transport configuration with a standard rib/spar wing- 
box. The layout of the rib/spar wingbox is shown in Figure |. Figure 2 presents the DLM aerodynamic panels used 
in this study. This structure and aerodynamic panel layout has been used in previous static aeroelastic optimization 
studies [7]. This model and a version with a higher aspect ratio, is being analyzed by the NASA Performance Adap- 
tive Aeroelastic Wing (PAAW) research project for a cruise condition at Mach 0.85. Dynamic aeroelastic optimization 
studies with the uCRM have been performed elsewhere, which included flutter optimization using DLM aerodynam- 
ics [24]. A similar optimization will be performed here, but with the URANS aerodynamics solved with the high- 
fidelity Navier-Stokes solver FUN3D v12.8. The Navier-Stokes code FUN3D (fully unstructured three-dimensional 
Navier-Stokes) is a finite volume unstructured CFD code for compressible flows [46, 47]. FUN3D solves the steady 
and unsteady Reynolds-averaged Navier-Stokes (URANS) flow equations loosely coupled with the Spalart-Allmaras 
turbulence model [48] on a tetrahedral mesh. The surface triangulation of this mesh is shown in Figure 3. System 
identification of the aerodynamic response to a perturbation of each mode is performed in the manner addressed in 
Ref. 33. The perturbation is a Gaussian pulse, which produces the generalized aerodynamic force response. A fast 


Fourier transform is applied to the time histories to produce the gAIC. 


A. Retained and Exanded Mode Sets 

Before discussing the flutter calculations, the influence of the number of modes in the expanded set in comparision 
to the retained number of modes will be investigated here. The cross correlation (Eq. 13) provides a metric to assess 
the adequacy of the expanded mode set chosen. The structural model has a complete set of 56242 nodes resulting 


in a large set of eigenvectors, whereas the expanded and retained mode sets are significantly smaller than this. Two 
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Fig. 1 uCRM wing-box structure from Ref. 7. 


Fig. 2 uCRM DLM panels from Ref. 24. 


sets of retained modes are constructed. One set has 20 retained modes (N,. = 20) and the other has 25 retained modes 
(N, = 25). The first set (NV, = 20) is constructed using no expansion of the mode set, i.e. N, = 20. The correlation 
measure is found in the second column of Table |. The correlation measure for the case in which N, = 25 while 
being constructed from a set with N, = 100 is found in the third column of Table 1. The remaining column gives the 
correlation measure for N,. = 20 and N, = 25. When N, = N, = 20, there is a marked drop off in the accuracy of the last 
8-9 modes. Retaining a much larger expanded set (NV, = 100) results in much better accuracy in all the modes, while 
in the case in which N, = 20 and N, = 25, there is some improvement in all the modes, and especially in modes 15-20. 


B. Flutter Onset 
Since the target cruise condition is Mach 0.85, attention to potential nonlinear transonic effects such as shock 
motion and shock separation is important to identifying flutter onset. Transonic nonlinearity can reduce the flutter onset 


Fig. 3 uCRM surface mesh for URANS calculation. 


dynamic pressure in a Mach range typically just below Mach 1. In the present configuration, nonlinear transonic effects 
are clearly evident at Mach 0.85. To see this effect, flutter simulations have been performed with DLM aerodynamics 
and URANS aerodynamics at Mach 0.70 and 0.85. The evolution of the modal flutter damping with increase in 
atmospheric density and speed are shown for the two Mach numbers in Figures 4 and_ 5. 

Using DLM aerodynamics, the onset of instability is generally around the same density at Mach 0.70 and 0.85. 
Using a URANS analysis, however, the onset of aeroelastic instability is clearly much lower at Mach 0.85. The 
significant difference between flutter onset using DLM and URANS aerodynamic models at Mach 0.85 motivates the 
present work. URANS flutter calculations using the method outlined in Sec. II will accurately provide flutter onset, 
while the approach of Sec. HI will be used to optimize the structure under the flutter onset constraint. 


C. Approximating the Flutter Onset of a Modified Structural Model Using Baseline DLM and URANS Aerodynamics 

To verify the performance of the method, the flutter onset of the baseline and modified finite element structural 
models will be compared. The gAIC matrix of the modified structure is calculated using the baseline and modified 
structural models and the URANS Q,. The aerodynamics are implemented into the flutter equations using a Roger 
approximation. Figure 6 shows the damping for the baseline and the modified models using URANS aerodynamics. 
These two models are sufficiently different to provide a realistic assessment of the current capability. 

The first step is to verify the model. The flutter onset for this configuration is a coalescing of the first three modes. 
Since these modes are the most significant, and differ most from the baseline to modified model, the gAIC matrix 
terms for the 1” two of these modes will be scrutinized. All the remaining results in this paper use N, = N, = 20. 
Figures 7 and 8 show four of the gAIC terms, each of which involves the second mode interacting with itself or one 
of the other first three modes. Figure 7 is derived from the DLM theory. Figure 8 is derived from URANS theory. For 
these gAIC terms, the baseline and modified models are easily differentiated. It is clear that the modified aerodynamic 
model using the method of Sec. II and that directly from the modified models are identical. This is true using either 
the DLM or the URANS theory. 


Table 1 £8 (Equation 13). 
mode N,=20 N,=25 N,= 20 
number N, = 20 N, = 100 N, = 25 


1 0.99999 1.00000 1.00000 
2 0.99833 0.99959 0.99931 
3 1.00023 1.00004 1.00010 
4 1.00034 1.00031 1.00030 
5 1.00042 0.99950 0.99981 
6 1.00090 0.99955 1.00020 
7 1.00112 1.00134 1.00130 
8 0.99985 0.99984 0.99999 
9 1.00092 0.99983 1.00110 
10 0.99860 1.00088 1.00130 
11 0.99011 1.00002 0.99881 
12 0.98912 0.99985 1.00320 
13. 0.99092 1.00041 0.99945 
14 0.98693 0.99883 1.00030 
15 0.97804 0.99890 0.99302 
16 0.98685 0.99814 0.99912 
17. 0.91360 0.99384 0.99718 
18 0.95756 0.99620 0.99481 
19 0.97184 0.99208 0.97882 
20 0.96091 1.00027 0.97577 
21 — 097525 — 

22 — 0.99766 — 

23 — 0.93356 — 

24 — 1.00524 — 

25 — 1.00055 — 
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Fig. 4 Flutter damping versus density for the baseline model using DLM and URANS aerodynamics, Mach 0.70, 
(O DLM, e URANS). 
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Fig. 5 Flutter damping versus density for the baseline model using DLM and URANS aerodynamics, Mach 0.85, 
(O DLM, e URANS). 
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Fig. 6 Flutter damping versus density using URANS aerodynamics, Mach 0.85, (e modified, © _ baseline). 


Figure 9 compares the evolution of the flutter damping for the URANS modified model using the method of Sec. 
II and that computed directly from a gAIC matrix calculated by a URANS system identification of the modified modes. 
Except for slight differences in some of the higher frequency modes, the two models are identical. The differences in 
the higher frequency modes are most likely due to a drop off in accuracy given that N. = N,. 

Finally, a comparision is shown in Figure 10 between the flutter damping calculated using the URANS baseline 
model, the URANS model created using the method of Sec. I, and damping calculated from a time-marching URANS 
flutter simulation. The damping of the time marching flutter solutions is obtained from the modal time histories by 


the log decrement method and flutter frequency by the spacing of zero crossings. The URANS time-marching flutter 
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Fig. 7 DLM gAIC terms, (— — — baseline, modified, @ modified using equations 8 and 11). 


simulation is started from a static aeroelastic solution at a slightly negative angle of attack in order to approximate 
a zero-lift condition. The URANS time-marching flutter simulation also includes wing weight while the URANS 
system identification solution did not include wing weight and is also started from an undeformed vehicle shape. 
Furthermore, the initial modal velocity used to start the time-marching flutter solutions was somewhat larger than 
the excitation used to create the analytical aerodynamic models. Because of these differences in the way damping 
and flutter frequencies are obtained, some differences are to be expected between the time-marching flutter and the 
frequency domain solutions. The time-marching solution damping (Re[p]) is calculated from the true damping (C) 
and the flutter frequency (@,) at two dynamic pressures near the flutter onset predicted by the present theory for both 
the baseline and the modified models. The damping calculated from the time-marching solutions is about a percent 
below the predicted flutter onset point in either case. On the other hand, the slope of the damping near flutter onset, 


and the increment in flutter onset from the baseline to modified solutions is well simulated by the present method. 


D. Optimization Results 

The optimization problem considered here seeks to find the structural sizing layout of the wingbox, which min- 
imizes the structural mass while still satisfying several constraints on aeroelastic performance. Design variables are 
based on the patch-level discretization in Fig. 1, where each skin panel is a design patch, each rib is a patch, and 
each spar section is a patch. Design variables include the shell thickness of each patch (bounded between 3 and 30 
mm), the thickness of the stiffeners attached to each patch (bounded between 2.5 and 30 mm), and the height of the 
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Fig. 8 URANS gAIC terms, (— — — baseline,— modified, @ modified using equations 8 and 11). 
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Fig. 9 Flutter damping versus density from URANS aero, Mach 0.85, generated with modified structure modes (e) and 
URANS transformed aero using equation 8 (C). 
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Fig. 10 Flutter damping versus density, comparing analytical with time-marching URANS. 


stiffeners attached to each patch (bounded between 30 and 150 mm). The stiffeners are modeled as smeared T-shaped 
stiffeners [49], where the flange is bonded to the shell members. The total number of design patches is 237, and the 
total number of design variables is 711. 

The objective function is the structural mass, computed by the volume of the wingbox finite element model. This 
function is minimized subject to many constraints. First, a flutter constraint is enforced at a Mach number of 0.85, 
across a range of flutter margins, using the methods outlined in Sec.s I and III. A series of static aeroelastic constraints 
are also enforced, spread across four trimmed maneuver loads: a 2.5g pull-up, a -1g push-over, a roll maneuver, and a 
2g landing load. All cases are run at the aircraft dive speed at sea level with full fuel. For each load case, the stresses in 
each finite element are computed, and the KS function is used to compress all of the elemental failure function values 
within a given patch into a single metric. 

Next, buckling analyses are run for each stiffened panel in the upper and lower skins using a Rayleigh-Ritz method 
(assumed buckling modes). Both global buckling of a stiffened panel (bordered by ribs and spars) and local buckling 
between each stiffener is computed, where simply-supported boundary conditions are used for both scenarios. As 
with the stress metrics, each buckling eigenvalue for a given stiffened panel is compressed into a single KS function. 
Design constraints are placed on both the stress and the buckling KS functions. The final gradient-based optimization 
problem is solved with the Globally-Convergent Method of Moving Asymptotes tool [50]. 

The optimal structural mass, as a function of the flutter margin constraint boundary (using a dive speed of 185 
m/s), is shown in Fig. 11, for two cases. The first case uses URANS to compute the baseline gAICs Q,, whereas 
the second case uses linear DLM aerodynamics. Both sets of results then use the transformation methods detailed 
above, to compute Qy as the structure is optimized. For low flutter margins in Fig. 11, the flutter constraint is very 
easy to satisfy, and is inactive (i.e., the static aeroelastic constraints drive the design process entirely). The structural 
mass (objective function) in this case is 8,631 kg. Increasing the flutter margin makes the constraint harder to satisfy, 
resulting in an increase in the structural mass. As also seen in the figure, the URANS model predicts a more aggressive 
flutter mechanism than the linear DLM model and therefore, results in heavier structural models. This added weight is 
due to the fact that the DLM does not accurately predict the drop in the flutter onset dynamic pressure typically seen 
in the transonic Mach range, frequently called the transonic dip. At a margin of 1.2, the URANS flutter model is 4.3% 


heavier. 
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Fig. 11 Optimal structural mass for various flutter constraint boundaries. 


The equivalent skin thickness distributions (sum of the panel thickness and the smeared stiffener thickness) are 
shown in Fig. 12, for three cases: the design with an inactive flutter constraint, and the two designs (URANS and 
DLM) with a 1.2 flutter margin constraint. For all cases, peak thickness tends to occur between the wing root (span 
location of 3 m) and the Yehudi break in the trailing edge (span location of 10 m), and then trends toward minimum 
gage at the wingtip. The addition of the flutter constraint increases the thicknesses in various areas, as expected, but the 
redistribution of material is quite different for the two flutter models. Particularly in the lower skin, the URANS model 
forces an increase in thickness localized near span location of 10 m, whereas the DLM model shows more of a global 
thickness increase across most of the skin. Such a result indicates that while the two aerodynamic models are both 
linearizations (in the case of DLM, the linearization is about the null state), and the flutter-damping plots presented 
above look topologically similar, the two flutter models in fact rely on inherently different flutter mechanisms, and 
therefore force the optimizer to react in fairly different ways. Neither model, it should be noted, is able to account for 
changes in the underlying shock structure with changes in the structural sizing variables, as the URANS linearization 


is conducted just once, and the resulting Q, matrix projected from one modal space to the next. 


V. Concluding Remarks 

A concise transformation of the unsteady generalized aerodynamic influence coefficient matrix based on the 
baseline structural modes is used to optimize a transport wingbox model under a transonic flutter constraint (among 
other linear static aeroelastic constraints), where the flutter constraint boundaries are computed using unsteady high- 
fidelity Reynolds-averaged Navier-Stokes aerodynamic model. Rather than recomputing a complete CFD-based time- 
marching flutter solution at each design iteration, the approach requires only a single matrix transformation. System 
identification techniques are used to extract linearized generalized aerodynamic influence coefficients for a baseline 
wing configuration. An inexpensive modal projection method is then used to transform these baseline aerodynamic 
terms into the modal space inhabited by the various structural iterates obtained during the optimization process. 

The results demonstrate the need for transonic flutter computations, versus purely-linear physics, and then demon- 


strate the accuracy and the efficiency of both the system identification techniques, and the modal projection techniques. 
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Fig. 12 Upper and lower skin equivalent thickness distributions. 


Finally, a structural optimization is presented with a large number of structural sizing design variables, which seeks to 
minimize structural mass under static aeroelastic stress and buckling constraints, as well as a flutter constraint. Opti- 
mization results are presented for a large range of required flutter margins, using both a CFD-based flutter constraint, 
and an entirely linear DLM-based flutter constraint. The CFD model predicts a more aggressive flutter mechanism 
than the linear DLM model (as the latter does not capture the flutter dip), and results in a heavier optimal structure. 
Furthermore, the optimal distribution of material under the flutter constraint identified by optimizing the structural 
models coupled with URANS and DLM aerodynamics is highly dissimilar. 
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